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Abstract 

A new model dissipation rate equation is proposed based on the dynamic equation of 
the mean-square vorticity fluctuation for large Reynolds number turbulence. The advan- 
tage of working with the vorticity fluctuation equation is that the physical meanings of the 
terms in this equation are more clear than those in the dissipation rate equation. Hence, 
the model development based on the vorticity fluctuation equation is more straightforward. 
The resulting form of the model equation is consistent with the spectral energy cascade 
an alysis introduced by Lumley (1992). The proposed model dissipation rate equation is 
numerically well behaved and can be applied in conjunction with any level of turbulence 
closure. In the present study, it is applied to a realizable eddy viscosity model. Flows that 
are examined include: (i) rotating homogeneous shear flows; (ii) free shear flows; (iii) a 
channel flow and flat plate boundary layers with and without pressure gradients; and (iv) 
backward facing step separated flows. In most cases, the present model predictions show 
considerable improvement over the standard k - e model. 
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1. Introduction 


The dissipation rate equation is often used in various turbulence closure schemes. Its 
exact equation can be written as, 

2u . 

e, t + Uie,i =ve.u - - —(p,kUi,k),i 

2l (1) 

- 2v 2 WJi 

where e = vvZJui~j and e' = vmjUij. All the terms on the right hand side of Eq. (1), 
except the viscous diffusion term ve.u, are new unknowns. Thus, they must be modeled 
before this equation can be used for applications. Modeling of these new unknowns, which 
are related the small scales of turbulence, is extremely difficult due to the lack of a clear 
understanding of their physical meanings. Therefore, in the literature, Eq. (1) has not been 
considered as a useful equation to work with. One often ignores Eq. (1) and creates a model 
dissipation rate equation which has a structure similar to that of the turbulent kinetic 
energy equation. In other words, the dissipation rate equation also has generation and 
destruction terms which are assumed to be proportional to the production and dissipation 
terms in the turbulent kinetic energy equation over the period of large eddy turn-over time, 
k/e. With this assumption, the resulting model dissipation rate equation can be written 
in the following form: 


e, t + Ui£ ti =ve t u - ( e'ui),i 


€ € 
—UiUjUij C £ 2~j^ 


( 2 ) 


Eq.(2) is the standard model dissipation rate equation which has been widely used in 
various turbulence closure schemes. In addition, several modified versions of Eq.(2) have 
also been also proposed for different applications, for example, in near-wall turbulent 
flows 1-4 and in rotating turbulent flows 5 . 


Recently, Lumley 6 proposed a dissipation rate equation based on the concept of non- 
equilibrium spectral energy transfer due to the interactions between eddies of different 
sizes. A new transport equation for an inverse time scale has also been suggested in 
conjunction with his new e equation which is of a different form from that of Eq.(2). This 
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model mimics the physics of the statistical energy transfer from large eddies to small eddies 
and is successful in the prediction of various turbulent free shear flows 6 . 

In this paper, we explore the possibility of deriving a new model form for the dissipa- 
tion equation which is not only more physically related to the original e equation but also 
simpler and more robust than the standard dissipation equation (2). This is achieved by 
first developing a model equation for the the dynamic equation for the mean-square vortic- 
ity fluctuation ZJl Wl. Once the dynamic equation for Wi Ul is modeled, a model dissipation 
rate equation can be readily obtained. 

2. Dynamic equation for W iui 


The exact equation for UiUi is 



),t + — v(— ),jj — 




-f" 

“f" "I" U1{U) jILi j 

VUi.jLOi.j 


( 3 ) 


where m and Ui are the fluctuating and mean velocities, and u ;» and Qi are the fluctuating 
and mean vorticities which are defined by 


LOi — £ijk U k,j i — £ijkUk,j ( 4 ) 

Tennekes and Lumley 7 clearly described the physical meaning of each term in Eq.(3). The 
first two terms on the right hand side represent the viscous transport and the turbulent 
transport of ZD~u57, respectively . The third term is the source term which is produced by 
fluctuating vortex stretching and mean vorticity. This term also appears in the equation 
for with the same sign, hence, it will either increase or decrease fijO; and UJUTi 

simultaneously. The fourth term represents the vorticity exchange between uJiuJi and 
because it appears with opposite sign in the equation for f UCti. The fifth term represents 
the source produced by mean vortex stretching. The sixth and seventh terms are the 
production due to fluctuating vortex stretching and the dissipation due to the viscosity 
of the fluid, respectively. Tennekes and Lumley 7 have shown that, at sufficiently higli 
turbulent Reynolds numbers, the sixth and the seventh terms are the largest terms in 
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Eq.(3). That is, 


zn/ W e>3/2' 1 

LOiUJjUij, PUi'jUJi'j ~ C'( j j 

The first, the third, the forth, and the fifth terms are the smallest, either of order (u 3 /! 3 ) 
or ( u/l) z Rf . In the above analysis, u O” denotes the order of magnitude. R t = ul/v is 
the turbulent Reynolds number, u and l are the characteristic velocity and length scales 
of turbulence, respectively. If the terms of order (u 3 jl 3 )Rt or larger were kept in equation 
(3), then the evolution of v&Ti should be described by the following equation, 

, UJiUJi V , uJTuJi, 1 f X ( 

( — 2 — )»* ~ + UiU} 3 u iJ (5) 

- vTdTjuTj 

As pointed out by Tennekes and Lumley 7 , at large Reynolds numbers, Eq.(5) becomes, 

Or equivalently, production equals dissipation. This relation indicates that the term 
UZJjUtj is always positive. In addition, it indicates that there is a new length scale 
created by the vortex stretching which is related to the derivative of fluctuating vorticities. 
The vortex stretching tends to reduce the size of eddies and to create a broad spectrum of 
eddy sizes. However, this process must end at a certain level of eddy size because of the 
smoothing effect of viscosity. We expect that the terminal eddy size is the Kolmogorov 
microscale which corresponds to the length scale of the derivative of fluctuating vorticity 
Uij. This can be easily verified from Eq.(6). 

In the following, a model equation for the mean-square vorticity fluctuation will first 
be developed. The proposed model dissipation is then derived straightforwardly. 


3. Modeling of the dynamic equation for UiUi 
3.1 Modeling of 

( 7 ) 

( 8 ) 


We first define a fluctuating anisotropic tensor using u>iUj 

COiUJ'i 1 - 

h u — _ J /)■ ■ 

0 — v 

J 1 . 1 , 1 . 1 , A 


VkWk 


then 


UiUjUij = bfjUkUkUij 
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We expect that the vortex stretching tends to align vortex lines with the strain rate and 
that the anisotropy b v. is mainly due to the fluctuating strain rate; hence, the anisotropy 
b?j may be assumed to be proportional to the strain rate Sij. That is, 

( 9 ) 


b?j oc 
3 s 


where 


This leads to 


S — (2 SijSijY^ i Sij — (UiJ ~i~‘Uj,i)/2 


SijUi r j , 


UJiUJjUij OC — = ^k^k S 

s 


If we further assume that UkWk and (^ s ij s ijY^ 2 are we ^ correlated, we may write 


Wi VjUij OC U k Uk -Jtoii'ii 

Noting that ZJiul = ZsipsTj at large Reynolds numbers, we may also write 


WjfcWfc UJiLOi 


(10) 


( 11 ) 


( 12 ) 


ViWjUi.j OC U) k Uk y/u)iWi = y-5 

J J y/UJiUfi 

Eqs.(ll) and (12) indicate that the model for uuuJuTj is of order (u z /£ z )R 3 t /2 as it should 


be. 


3.2 Modeling of UiUjUij - WijWij 


Eq.(2) indicates that - vJSi jUTj must b e of or der (u z /t z )R t . Therefore, 

—va}i_jU)ij must cancel oo^k (or uJ^jJk \J^sijSij) m suc ^ a way that their 

j i t*y 

difference is smaller than UJl uJZ UiuJi/ y/u>iu>i (or UkUJk y^SijSij ) by an order of R t . This 
suggests that the sum of these two terms can be related to the following two terms: 


UkUk S , 




(13) 


since both the ratio of kjv to \/uiLUi and the ratio of s to S are of order R 1 / 2 ■ k(~ u 2 ) 
denotes the turbulent kinetic energy and S is the mean strain rate (y/2SijSij). As a result, 
the dynamic equation for fluctuating vorticity can be modeled as 

(^r ),< + + Ciwxs 

2 2 1 (14) 

U>k<*>k UjUj v ' 

£ + y/uiidi 
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Here, the term ZJk&kS represents the part of the production term Z^uJuTj due to the mean 
strain rate which is always positive if 5 0. Note that the sum of last two terms in Eq.(14) 

models the last two terms in Eq.(6) as a whole and should not be viewed as a model for 
either individual term. 

4. Modeling of the dissipation rate equation 

Noting that at large Reynolds number ZJjuTi = uijiilj and e = uuiUi, we readily 
obtain a modeled dissipation rate equation, 

+ Uj cj = + Cl s S - c, - (15) 

where the model coefficients, C\ and C 2 , are expected to be independent of the Reynolds 
number as the Reynolds number becomes large. We note that C\ and C 2 may be affected 
by the solid body rotation imposed on turbulence through the reduction of fluctuation 
vortex stretching, uJiuJJufj, as was shown by Bardina®; however, this effect is rather weak 
compared to the other mechanisms. For example, Reynolds stresses will be first substan- 
tially affected by rotation and result in a substantial change of the turbulent field, say k, 
as shown in the calculation of the rotating homogeneous shear flows in section 5.1. This 
will also affect the evolution of e through, say, k. The signs of C\ and C 2 can be easily 
determined. For example, in a decaying grid turbulence, only the last term on the right 
hand side of Eq.(15) is non-zero and must be negative, hence C% must be positive. For the 
case of homogeneous shear flow, both the turbulent kinetic energy and its dissipation rate 
increase with time so that the source term in Eq.(15) must be positive hence 0\ must be 
positive. In fact, these two types of flows®’® will be used for calibrating the values of CJ\ 
and C 2 . 

The differences between the present model dissipation rate equation, Eq.(15), and the 
standard model dissipation rate equation, Eq.(2), are the following. First, the destruc- 
tion” term is well behaved so that Eq.(15) does not have a singularity anywhere in the 
flow field. This is due to the presence of -Jve. Although y/us is a term of lower order 
than that of k at large Reynolds numbers, it is kept here to avoid the singularity in case 
k — ► 0. Second, the Reynolds stresses do not appear in Eq.(15). Consequently, the present 
model dissipation rate equation will be more robust than the standard model dissipation 
rate equation when it is used in conjunction with second-order closure schemes, since S 
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normally behaves better than the Reynolds stresses in numerical calculations, especially for 
cases with poor initial conditions. In addition, the present form of the ‘ production term 
is similar to that proposed by Lumley 2 which is based on the concept of spectral energy 
transfer. It is believed that the present model describes the turbulent vortex stretching 
and dissipation terms more appropriately. 


Eq.(15) can be applied in conjunction with any level of turbulence closures. However, 
the turbulent transport term (Fui) t i needs to be modeled differently at different levels or 
turbulence modeling. Here, we apply Eq.(15) to a realizable eddy viscosity model of Shih 
et al. 10 , and is modeled as 





(16) 


For completeness, the realizable eddy viscosity model will be briefly described in the fol- 
lowing section. The model coefficients C\ , C 2 and will be determined later. 


4.1 A realizable eddy viscosity model 

Shih et al. 10 proposed a realizable Reynolds stress algebraic equation model. Here, 
for simplicity, we adopt its linear form: 


where, 


and 


—UiUj = vt (Uij + Uj t i) — —kSij 

_ r *1 r 2 / 3 

Vt C,l e' ,l 1.25 + 77 + 0.9f 


v = — , « = — . s = ^ 2 S^sr it n = Jw/n, 
£ ^ 


kCl 

e 


Sii — — (Uij + Ujj), Cl*j — — (UiJ Uj.i) + 4 ernjtWrn 

where o>„. is the rotation rate of the coordinate frame. 


(17) 

(18) 

(19) 

( 20 ) 


The modeled k and e equations are 

, vt 


k,t + Ujkj — (— kj\j UiUjUij e 

e.t + Ujej = (^j) J + Ci S c - C, k + ^ 


( 21 ) 

( 22 ) 
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where Ci, C 2 , &k and a e axe model constants. 


4.2 Model constants C\ , C 2 and a E 


In decaying grid turbulence at large Reynolds number, the equations for turbulent 
kinetic energy k and its dissipation rate s are 



£ ,t ~ 



(23) 


K / l \ — n c / \ — « 

5 - v ’ £0 ~ <0 

The following equations can be obtained from the k and e equations: 

a = n + l ,C 2 =^±i (24) 

n 

Experiments 8 show that the decay exponent n varies from 1.08 to 1.30. In this study 
we choose C 2 = 1.9 which corresponds to n = 1.11. After C 2 is chosen, we use the 
experimental data for homogeneous shear flow® to calibrate the value of Ci which is found 
to be about 0.42. 


The value of a e will be estimated by using the log-law in turbulent boundary layers. 


These relations are: 


lUIiog^ + c 


U-r 


K 


V 


—UV ~ UZ., 


— UV 


_dU 

dy 


which lead to the following, 


k 

u 


_ 1 dU _ u T 

l ~ y/c? d y ~ K y 

, kdu k 1 

V ~ S — „ a„. — „,2 


£ = 


Ky 


e dy u\ y[Ul 
Substituting the 77 and £ above into the Eq.(18) for C fl , we obtain 


(25) 


(26) 


1.25 C, t + 1.9 yc; - - = 0 


(27) 
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which gives C, t = 0.0864 in the log-law region. This value is then used below to determine 

o e . 

Analyzing the dissipation rate equation in the log-law region, we obtain 

/c 2 

ae ~c 2 ^c;-c l (28) 

= 1.21 

where the von Karman constant k = 0.41. The model constants are summarized in the 
Table 1. 


Table 1. Model constants 


Ok 

O e 


c 2 

1.0 

1.21 

0.42 

1.9 


5. Model applications 

The results of some turbulent flow calculations using the proposed new turbulence 
model, Eqs.(17)-(22), are shown in this section. These include (i) rotating homogeneous 
shear flows, (ii) free shear flows, (iii) channel flow and boundary layers with and without 
pressure gradients, and (v) backward-facing step flows. The computations with the present 
and the standard k - e models are compared with DNS, LES and various experiments. 

5.1 Rotating homogeneous shear flows 

The comparisons are made with the large eddy simulation of Bardina et al . 5 for foui 
different cases of Cl/S (which are Cl/S= 0.0, Cl/ S=~0.50, Cl/ S— 0.25, and Cl/S 0.50). The 
initial conditions in all these cases correspond to isotropic turbulence and eo/Sk 0 = 0.297. 
Figure 1 (a) compares the evolution of turbulence kinetic energy, normalized by its initial 
value fc 0 , with the non-dimensionalized time St for Cl/S = 0.0. For this case both the 
present and the standard k — e (denoted by ske hereafter) models show the trends exhibited 
by LES, with the present model being closer to the LES data. Figure 1 (b) shows the 
comparisons for the case Cl/S = 0.25. The ske model does not pick the effect of rotation 
on turbulence as it gives the same growth rate of turbulence kinetic energy as it did for the 
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no rotation case, a result which is already known. On the other hand the present model is 
in reasonable agreement with the LES data as it shows a faster growth of the turbulence 
kinetic energy with time over the no rotation case. Figures 1 (c) and 1 (d) compare the 
evolution of turbulence kinetic energy for two more cases of fi/S = 0.5 and £l/S = -0.5. 
For the first of these cases the LES shows that the growth rate of the turbulence kinetic 
energy is decreased over the no rotation rate case. The present model is able to pick up 
this trend and although the agreement between the present model and the LES is not as 
good as it is for the other cases, it still is a lot better than the ske model. For the case or 
0/5 = -0.5 the agreement between the present model and the LES is quite reasonable. 
It should be pointed out that the rotation effect on turbulence is brought about by the C, t 
formulation of Eq. (18). 

5.2 Free shear flows 

Calculations using the present and the ske models were performed for a mixing layer, a 
plane jet and a round jet. Figures 2, 3 and 4 show the comparisons of the self-similar profiles 
between the model predictions and the measurements for the mixing layer, plane and round 
jets, respectively. In these figures, the profiles for the mean velocity the Reynolds shear 
stress and the turbulent kinetic energy are presented. For the mixing layer, the results are 
shown in a self-similar coordinate rj defined as 

y — Jfo.5 

V • 

2/0.9 - y o.i 

Vo.i , 2/0.5. and yo.g denote the locations where the ratio of the local mean velocity to that of 
the free stream are 0.1, 0.5, and 0.9, respectively. Figure 2 shows that the mean velocity 
profiles of the mixing layer predicted by either the present model or the ske model agree well 
with experimental data of Patel 11 . The present model, however, gives better predictions 
of the turbulent kinetic energy and the Reynolds shear stress distributions than the ske 
model. This is especially true for their peak levels. The predictions for the plane jet 
are shown in Figure 3. The model predictions are compared with the measurements of 
Gutmark and Wygnanski 12 , Bradbury 13 , and Hekestad 14 . The predictions given by both 
the present model and the ske model agree well with the experimental data. The turbulent 
kinetic energy level at the jet centerline is slightly lower than the measured values. For the 
round jet, the comparisons are made between the model predictions and the measurements 
of Wygnanski and Fielder 15 and Rodi 16 and are shown in Figure 4. The profile distribution 
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of the mean velocity predicted by the current model agree well with the Rodi’s data, while 
the ske model predicts a faster spreading of the round jet into the surroundings and a wider 
distribution. Significant improvement is also achieved in the prediction of the turbulent 
kinetic energy profile over the ske model in terms of both the centerline level and the 
overall distribution. The new modeling of the production of the dissipation rate according 
to its dynamic significance is believed to contribute to these improvements over the ske 
model. The calculated spreading rates of these flows are compared with measurements 
and are shown in Table 2. 

Table 2. Comparison of the spreading rates of turbulent free shear flows 


Case 

measurement 

ske model 

present model 

mixing layer 

0.13-0.17 

0.152 

0.144 

plane jet 

0.105-0.11 

0.109 

0.103 

round jet 

0.085-0.095 

0.116 

0.105 


The present model demonstrates better predictions over the ske model; especially, the 
spreading rate anomaly of plane and round jets are significantly reduced. This improvement 
is mainly due to the new dissipation rate equation. 

5.3 Channel flow and boundary layer flows 

Turbulent channel flow and turbulent boundary layer flows with/without pressure 
gradients were calculated to test the performance of the present model for wall bounded 
flows. Since the present model is proposed for turbulent flows away from the wall, the 
integration was carried out down to = 80 rather than to the wall in the calculations. 
At y + = 80, DNS values were used as the boundary conditions for the turbulent channel 
flows and wall functions were used for the turbulent boundary layer flows. 

The results for 2D fully developed turbulent channel flow at Re r = 395 are shown in 
Figure 5. This flow was calculated by Kim 17 using direct numerical simulation. Both the 
present model and the ske model agree reasonably well with the DNS data. Figure 6 shows 
the results for the flat plate boundary layers with the Reynolds number up to Reo = 16000. 
Here, comparison is made with the experimental results of Wieghardt 18 . Both the present 
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model and the ske model give good agreement with the experiments. Overall, the present 
model gives a slightly better prediction for boundary layer development. 

Figures 7 shows the results for the Herring and Norbury flow 19 , which is the bound- 
ary layer flow under favorable pressure gradient. The present model and the ske model 
give similar predictions for the velocity profile. The ske model gives better predictions for 
the skin friction coefficient and the shape factor while the present model gives a better 
prediction for the boundary layer thickness. The turbulent boundary layer under adverse 
pressure gradient studied by Bradshaw 29 and the turbulent boundary layer under increas- 
ingly adverse pressure gradient studied by Samuel and Joubert 21 were also calculated. The 
results are shown in Figure 8 and Figure 9, respectively. In both cases, the present model 
gives better predictions, particularly for the skin friction. 

5,4 Backward-facing step flows 

The performance of the present model for complicated recirculating flows is demon- 
strated through calculations for two backward-facing step flows, one (DS-case 22 ) with 
smaller and the other (KKJ-case 23 ) with larger step height, both of which have been 
extensively used to benchmark calculations of separated flows. The calculations were per- 
formed with a conservative finite- volume procedure. The convection terms of the governing- 
equations were discretized by a second-order accurate and bounded differencing scheme , 
and all the other terms by the standard central differencing scheme. Sufficiently fine grids, 
with 201x109 points in the DS-case and 199x91 points in the KKJ-case, were used to 
establish numerical credibility of the solutions. The computational domain had a length 
of 50 step heights, one fifth of which was placed upstream of the step. The experimental 
data were used to specify the inflow conditions and the fully-developed flow conditions 
were imposed at the outflow boundary. The standard wall function approach 5 was used 
to bridge the viscous sublayer near the wall. The comparison of the predicted and the 
measured reattachment length is shown in Table 3. Figures 10-14 show the comparsion of 
the skin friction and pressure distributions along the bottom wall, and the mean velocity 
as well as the turbulent stress profiles at three downstream locations. All the quantities 
were normalized by the step height h and the experimental reference free-stream velocity 

Urr.f • 
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Table 3. Comparison of the reattachment point locations 


Case 

measurement 

ske model 

present model 

DS 

6.1 

4.99 

5.36 

KKJ 

7± 0.5 

6.35 

6.8 


It can be seen that the overall performance of the present model is better than that 
of the ske model. 


6. Concluding Remarks 

A new dissipation rate equation based on an analysis of the dynamic equation for 
fluctuating vorticity is proposed. The form of the model equation is consistent with the 
one based on the spectral energy cascade model introduced by Lumley 6 . The proposed 
dissipation rate equation can be readily applied to any level of turbulence closures. In the 
present study, it has been applied to the new eddy viscosity model described in Section 
4.1 and tested in various benchmark flows including: rotating homogeneous shear flows; 
boundary-free shear flows; channel and flat boundary layer flows with and without pres- 
sure gradients; and backward facing step flows. The results show that the present model 
performs better than the standard k-e model in almost all of the cases tested. The new 
model dissipation rate equation is expected to enhance the numerical stability of turbu- 
lence model equations, especially when it is used in conjunction with the advanced closure 
schemes, such as second order closures. We have just finished implementing the present 
model into the LRR 26 second order closure. Preliminary results seem to show that the 
initial decay behavior of and the effect of rotation on both k and e of initially isotropic 
rotating homogeneous shear flows are well captured. 
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Figure 2. Self-similar profiles for a plane mixing layer, (a) mean velocity; (b) turbulent kinetic energy; (c) 
Reynolds shear stress. 
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Figure 3. Self-similar profiles for a plane jet. (a) mean velocity; (b) turbulent kinetic energy; (c) Reynolds shear 
stress. 
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Figure 4. Self-similar profiles for a round jet. (a) mean velocity; (b) turbulent kinetic energy; (c) Reynolds 
shear stress. 
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Figure 6. Zero pressure gradient turbulent boundary layer, (a) mean velocity at Re$= 8900; (b) skin friction 
coefficient; (d) displacement thickness; (d) momentum thickness. 
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Figure 11. Static pressure coefficient C p along the bottom wall 
(legend as in figure 10) 
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Figure 12. Streamwise mean velocity Z7-profiles 
(legend as in figure 10) 
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Figure 13. Turbulent shear stress profiles in the DS-case 
(legend as in figure 10) 
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Figure 14. Turbulent normal stress profiles in the DS-case 
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